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Abstract 

Dendrites are extensions to the neuronal cell body in the brain which are posited in several functions 
ranging from electrical and chemical compartmentalization to coincident detection. Dendrites vary across 
cell types but one common feature they share is a branched structure. The cable equation is a partial 
differential equation that describes the evolution of voltage in the dendrite. A solution to this equation 
is normally found using finite difference schemes. Spectral methods have also been used to solve this 
equation with better accuracy. Here we report the solution to the cable equation using a compact finite 
difference scheme which gives spectral like resolution and can be more easily used with modifications to 
the cable equation like nonlinearity, branching and other morphological transforms. Widely used in the 
study of turbulent flow and wave propagation, this is the first time it is being used to study conduction in 
the brain. Here we discuss its usage in a passive, soma dendrite construct. The superior resolving power 
of this scheme compared to the central difference scheme becomes apparent with increasing complexity 
of the model. 
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Introduction 

The cable equation governs the dynamics of membrane potential in thin and long cylinders such as axons 
or dendrites in neurons. It is Wilfred Rail's derivation of the cable equation which set the ground for 
describing quantitatively the change in voltage across space and time in dendrites [I]-[3]. Under consid- 
eration here is a linear system involving neuronal processes that have a voltage independent component. 
The model is a passive one which means that membrane conductance is fixed. Here we only look at 
subthreshold responses which are about half the strength required to generate action potentials. If the 
response is stronger the membrane potential will begin to change and our assumption of passive dendrite 
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will no longer be true. The cable equation is ( [2J, equation 2.7) 

- (V m (x,t) - Vrest) + r m Iinj(x,t) (1) 



r m dV m (x,t) _ X 2 d 2 V m (x,t) 



at ax 2 

Here V m is the membrane potential in mV, r m is the membrane resistance per unit length of the fibre 
in fi.cm. Ii n j(x,t) is the injected current in Amperes, r m = r m c m is the membrane time constant and 
A = (fm/fa) 1 is the membrane space constant. c m is capacitance per unit length of cable of diameter d 
in units of F/cm. r a is axial resistance per unit length in J7/cm. V re st varies between —50 and —90 mV 
depending on cell type. Nondimensionalising this equation results in : 

dV(X,T) _ d 2 V{X,T) I(X,T) 

or - ox 2 ~ v{x ' T)+ ^r (2) 

Here X = x/X, T = tr, I{X,T) = \Tl in j(x,t), I in j(x,t) is stimulus current density and c m = C m ird, C m 
is specific capacitance/unit area in F/ cm 2 . Fig [T] shows the compartmcntal electrical representation of 
a segment of passive membrane. 



Numerical method for the cable equation 

Numerical analysis involves designing algorithms that approximate the solution to the equation under 
consideration As in all approximations, this results in a small error called the truncation error. Most 
finite difference schemes use the central difference formulae. It can be derived from the Taylor's series 
and shows algebraic convergence [5]. The second derivative ^£ in equation 1, also written as V" can be 
approximated as ( [5], equation 8.54) : 

VI' = W+i- 2 ^ + ^-i) + {h% 2 < i < N 1 (3) 

The subscripts denote locations on a uniformly spaced grid Xi and h is the mesh size Xi+i — Xi. 0{h 2 ) is 
the truncation error which is second order in this case. The length L of the cable is divided into uniform 
segments of length h = L/(N — 1) where N is the number of nodes. Vi is the voltage at location i at 
distance Xi = (i — l)h. Indices i = 1 and i = N denote boundary points X = and X = L. Fig [2] shows 
this. 
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For spectral methods the convergence is exponentially fast [BJ. The cable equation has 



been solved using a spectral Chebyshev method [7] . The solutions have uniform, high numerical accuracy 
at any spatial point- not just the original collocation points. The truncation error is negligible and the 
total error can be the roundoff error of the computations. 

In this paper a compact finite difference scheme is utilized [8]. Such a scheme provides better resolution 
at the shorter length scales compared to the usual explicit schemes like equation 3. Thus it provides the 
resolution of spectral methods but unlike spectral methods which are sensitive to discontinuities, requires 
little extra procedure for generalisation of the governing equation ( like geometry and other property 
variations along the dendrite). 

Spatial discretisation : Using compact finite difference schemes to solve the cable equation 

The second derivative is approximated using the following equation ( [8], equation 2.2): 



where V" represents the finite difference approximation to the second derivative at node i. The relations 
between the coefficients a,b,c and a , 8 are derived by matching the Taylor series coefficients of various 
orders. We take ( [SJ, equation 2.2.7) 



BV('_ 2 + aVjU + V{' + aV/li + 8V, 




c{V l+3 - 2V + Vj- 3 ) 
9h 2 

b(V i+ 2 - 2Vj + Vj_ 2 ) 
4h 2 

a(V i+1 -2V t + V^ 1 ) 
h 2 



,2<i<N-l 



(4) 



a 



— ,8 = 0,a = — ,b = — ,c = 
IV H 11 11' 



to obtain a sixth order formula. For the boundaries the scheme chosen is ( [8], equation 4.3.4) 



V" + aVz 



aVi + bV 2 + cV 3 + dVj + eV 5 
h 2 



(5) 
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A similar equation connects and V^_ 1 . By requiring third - order formal accuracy the coefficients 
are reduced to ( jS], equation 4. 3. 6), choosing a = ^ from classical Pade scheme which is fourth order. 

(11a + 35) , -(5a + 26) (a + 19) , (a - 14) 

a = , o = , c = , a = 

12 3 2 3 

(11 - a) 
e= ^2— 

The truncation error is reduced to ( a ~ 2 10 )fe 3 V^ 5 . If a is 10, truncation error becomes h 4 . Equation 4 and 
5 applied at interior points result in a matrix problem AV" = B where A is tridiagonal and V" can be 
obtained easily. 



Time discretisation 

The values for V" calculated from the compact-difference scheme were used to integrate the result in time 
using an explicit time stepping scheme - forward Euler. If t = nAT, V n = V(t) and V n+1 = V(t + AT), 
then : 

yn+l = yn + y.£yn nAT)AT (6) 

Stability conditions have dictated the choice of the time step 

Ax 2 

AT = ^- (7) 

It varies as shown in (Table [1]). The numerical integration in time has been done with an explicit scheme. 
Since spatial derivatives are obtained with a compact scheme, which is an implicit formula that requires 
the solution of a linear system, implicit time-stepping is not possible. Implicit time-stepping is desirable 
to overcome the severe restrictions that stability imposes on the time-step of conditionally stable explicit 
schemes. A work-around is to use a predictor-corrector scheme which uses an explicit step estimate from 
the predictor step in a corrector step which is also an explicit step. 

Computations were performed on a Toshiba Satellite Pro laptop using Octave in a Linux(Ubuntu)cnvironment. 
The data used for simulations is given in (Table [3]) . 
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Specific problems 

There are four problems we have considered. Before discussing each in detail, the initial and boundary 
conditions used are discussed: 



Initial Condition: 



V{x,0) = v(x),0 < x < I (8) 



Here v(x) — — 70 mV. 
Boundary Condition: 
Notation : V'(x,t) = 
Current Injection: 

If there is a current of magnitude I(t) injected into the end at x = ( in the positive x- direc- 
tion)^ [3J, equation 4.36) 

V'(0,t) = -nl(t),t>0 (9) 
If the current I(t) is injected into the end at x = l,( [3], equation 4.37) 

V'(l,t)= ri I(t),t>0 (10) 



Sealed end 

If end at x = is sealed, ( [3], equation 4.32) 

V'(0,t) =0,t > (11) 
If end at x = I is sealed, ( [5], equation 4.33) : 

V'{l,t) = 0,t > (12) 

Killed end 

If end at x — is killed, ( [3], equation 4.34) 



1/(0, t) = 0,<>0 



(13) 
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or if end at x = I is killed, ( [3J, equation 4.35) 

V(l,t) = 0,t>0 (14) 

1. Point soma dendrite, current injection at x = 0: The point soma dendrite construct we have 
used is a soma( cell body) collapsed into a point as mentioned in [7] along with a dendrite attached to 
it (Fig. U). The governing equation is the cable equation as mentioned above in equation 1. ( Range for 
independent variables < x < l,t > and the steady state at t — > oo). Initial condition is given by 
equation 8. Boundary conditions are current injection at i = 1 (equation 9) and at i = N, sealed end 
(equation 12) or killed end(cquation 14). 

The analytical equation at steady state for this formulation with sealed end boundary condition at i = N 
is ( [3J, Table 4.2) : 

where V(x) is the voltage at any given i, V(0) is the voltage at i = 1 , L = I, X = The analytical 
equation at steady state for the boundary condition with killed end at i = N is ( [3j, Table 4.2): 

n*) = v m Sin £%- L p ) (i6) 

where V(x), V(0),L and X are as defined above. 

2. Lumped soma model, current injected into soma: Additionally, we have used a lumped soma 
model where the current is injected into the soma (Fig. 2]). Here the boundary condition at i = 1 is a 
lumped soma boundary condition at steady state ( [3J, equation 6.14), For a first order case, this is : 

=V(o) - 7 V'(0) (17) 

T 

where Rs is the membrane resistance at soma, 7 = I s is the current injected at the soma, r is the 
time constant. 

The boundary condition at i = TV can either be sealed end ( equation 12) or killed end ( equation 14). 
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The analytical equation at steady state is ( [3J, equation 6.20) : 

where V{X) is the voltage at any given node i, I is the injected current, L = j, X ~ j and = -^- + 7^7 
( [3J, equation 6.21). 

3. Point soma dendrite, current injection at i = N: When current is injected at the end of the 

dendrite at i = N as shown in (Fig. 2]) , then according to convention longitudinal currents are positive in 
the positive x- direction. The equation governing this is the cable equation (1) as before with conditions 
given earlier. The boundary condition for current injected at the end (i = N) is defined by equation 10. 
The boundary conditions at the end (i = 1) can be sealed ( equation 11) or killed ( equation 13). 

4. Point soma dendrite, current injection at i = N/2: Here current is injected at the point i = N/2 
as shown in (Fig. 

The current injection condition is 



I(X,T) = \ T I inj {x,t),i = N/2 (19) 

I(X,T) = 0,i ^ N/2 (20) 

The boundary conditions at i = 1 can be sealed ( equation 11) or killed (equation 13). Similarly the 
boundary conditions at i = N can be sealed ( equation 12) or killed ( equation 14). 
The aim of this paper is to solve the cable equation using the compact difference scheme. This has been 
achieved in the passive soma dendrite construct. To the best of our knowledge, this is the first time the 
compact difference scheme is being used to solve equations involving changing parameters in the brain. 



Results 



The results obtained in the above discussed simulations are presented here: 
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Point soma dendrite, current injection at x = 
Current injection at i = 1, Evolution of voltage 

In this simulation, current is injected at i = 1 and the evolution of voltage at various times is observed. 

i = (AT)(n)(r) (21) 

For N = 30, r = 20 msec,the voltage evolution is looked at t as shown in (Table [2]). The results for both 
sealed and killed end boundary condition at i = N can be seen in (Fig [5j. 

Current injection at i = 1, Comparison with analytical equation at steady state 

In this simulation, current is injected at i = 1, N = 30 and the simulation is run till t = 500 msec. The 
results are compared with the analytical equation. V/Vq is plotted against X. The results are shown in 
(Fig.®. 

Current injection at i = 1, Comparison between compact, central difference schemes with 
analytical equation at steady state, sealed end boundary condition 

y- is plotted against X for the compact difference scheme, central difference scheme and the analytical 
solution for sealed end and results are shown in (Fig. [7]) . The second derivative is approximated 
using equation 4. By choosing coefficients as given below, it gives fourth order accuracy. The coefficients 
are ( [3], equation 2.2.6): 

r n n 4 (!-") h l(-l + 10tt) 1 , 1 

/3 = 0, c = 0, a=^— , b= , a=-, al = - 

The equation for the fourth order central difference scheme is calculated using ( [5], Table 6.4): 

V"(i) = ^-- ■ 2—^ ^—i- 2— ^— -,2<i<N -1 (22) 

The boundaries for the central difference scheme are calculated using : 

V"(l) = ( 4V(l + l)-3y(l)-V(l + 2) ) (23) 
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v» { N) = ^-D-sy-m-^ (24) 

The points next to the boundary are calculated using : 

v , l{2) = fX(g) - - 2) + ^F(4) - 13VX5) + §V(6) - f V(7) (2g) 

^ _ f y(7V - 1) - fV^V - 2) + ^y(7V - 3) - 13V (N - 4) + §^(7V - 5) - |VXjV - 6) 



(26) 



where h = Ax. 

Comparison of error between compact and central difference schemes, N = 10,20,30,40 

%Errorl = 100 * | Vcomp ~^ a " al | (27) 

%Error2 = 100 * | ^% ~ K " Qi | (28) 

Vana/(1) 

%Error\ and %Error2 are plotted against X. The results are shown in (Fig. |8]). 
Dependence of voltage changes on dendritic diameter 

Current is injected at i = 1 for dendrites of various diameter and y- is plotted versus X. The results are 
given in (Fig. [HJ) - 

Point soma dendrite, current injection as x = I 

Current injection at i = N, sealed and killed end, Evolution of voltage 

The current is injected at the end opposite to the soma and the voltage evolution from t = 0.0016969msec 
to t = 500msec is plotted. The intial condition is shown in equation 8. The boundary condition at i = 1 is 
chosen to be either sealed end given by equation 11 or killed end, equation 13. The boundary condition at 
i = N is the current injection boundary condition shown by equation 10. Results are shown in (Fig. [TT 



Point soma dendrite, current injection as x 



i 



The current is injected at i = The initial condition is given by equation 8. The boundary conditions 
at i = 1 and i = N can either be sealed end given by equations 11, 12 or killed end given by equations 
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13,14. 

Lumped soma, current injection at soma 

Current injection at soma, lumped soma boundary condition, Evolution of voltage 

Here the current is injected at the soma as shown in (Fig. 0]) and the evolution of voltage at various times 
is observed. 

t=(AT)(n)(r) (29) 

For N = 30, t = 20 msec, the voltage evolution is looked at t as shown in (Table [2]). The results for both 
sealed and killed end boundary conditions at i = N can be seen in (Fig[5j). 

Discussion 

It might be apt here to quote Wilfred Rail [TO],' An important basic principle of prudent research, both 
theoretical and experimental, is not to tackle too many complications at once. This was my reason for 
beginning my modeling with uniform, passive membrane, and with idealized geometry. Once you solve 
the reduced problem, you can then begin to deal with some of the complications judged to be functionally 
important. ' Keeping this advice in mind, we have demonstrated here the usage of the compact difference 
scheme to a passive soma dendrite construct under varying experimental and morphological situations. 

Our results show that the scheme is robust under all the conditions and shows a good fit 
with the analytical solution even at N = 10. By looking at the percent error EP1 ( point soma), (Table. 2|) 
at N = 10, it can be seen that it is almost the same for compact and central schemes. The error due to 
the central scheme is about 0.001 less than that due to the compact scheme. 

Lele [5] defines resolving efficiency as a difference between the modified wavenumbers of 
spectral scheme and the differencing scheme. This is specific to a scheme for any given error. Thus it 
can be seen that at error e <= 0.001 , the resolving efficiency of the fourth order compact scheme is 0.22 
which is greater than that of the fourth order central scheme which is 0.17. Sec (Table. [S]). It is also seen 
that the resolving efficiency of the sixth order compact scheme is 0.38. 

The simulations dealing with evolution of voltage show that when current is injected at 
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i = l, the voltage evolves to a higher value for both sealed and killed end. The cable reaches a steady 
state faster for the killed end than the sealed end. This is to be expected as there is no leakage at the 
sealed boundary. In the lumped soma case, the voltage reached at steady state is less than that at the 
point soma. This could be due to the R-C circuit comprising the lumped soma which results in some 
current injected going to charge this resulting in a smaller value of current flowing through the cable. 
The shape of the voltage evolution when current is injected at i = Appoint soma) is different for both 
sealed and killed end and it also evolves to a lesser voltage at steady state than current injection at i = 1. 
The rate of evolution is also much slower in this case. Finally, when current is injected at i = the 
resulting plots show predictable curves expected when either both ends are sealed or killed. 

As A is directly proportional to the square root of the diameter of the dendrite, one can 
see that as the diameter increases (Fig [9J , the lambda increases too and the plots for both sealed and 
killed end boundaries give expected results. See (Table. 

The numerical integration under consideration is an explicit scheme. In compact schemes, 
explicit schemes can only be used. This is unlike central schemes where both explicit and implicit 
techniques can be used. In the compact scheme used here, the stability is linked to AT or the time step 
that is being used. Various values of AT used are shown in (Table [1} . It is observed that this results in 
slower computations. 

In the simple cable chosen here, it is not possible to see the spatial resolving efficiency of 
the compact scheme. It is expected that the spatial resolving power of the compact scheme can be viewed 
with greater complexity of the model and also by choosing schemes with different coefheiants than the 
sixth order tridiagonal scheme chosen here. From Lele [5] it can be seen that an eighth order tridiagonal 
scheme can yield at e <= 0.001, a resolving efficiency of 0.48 and a spectral like pentadiagonal scheme 
can yield a resolving efficiency of 0.84. 

In this paper, it has been shown that the compact scheme can be used to solve the cable 
equation under various morphological and experimental conditions. Ongoing work looks at the use of 
this tool in the case of branched, tapered and active dendrites where the complex geometry could yield 
situations where the compact scheme is more efficient. 
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Figure 1. Compartmental electrical representation of a segment of passive cable 



- sealed/ killed end 



node 
pointsoma 



dendrite 



Figure 2. Discretisation of a dendritic tree 
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Figure 3. Dendrite with lumped soma 
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-Currentlnjection ' 



Dendrite 
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(a) Point soma dendrite construct (b) Soma dendrite construct with 
with current injection at i = 1 lumped soma boundary condition 

Currentlnjection , Currentlnjection 



(c) Point soma dendrite construct (d) Point soma dendrite construct 
with current injected at i = N with current injected at i = N/2 



Figure 4. Soma dendrite construct - different cases 
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(c) sealed end, lumped soma (d) killed end, lumped soma 



Figure 5. Evolution of voltage point and lumped soma, current injected at i = 1, 

n=l; — : n = 100; :n == 500; :n == 2000; :n == 294670 
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(a) sealed end (b) killed end 



Figure 6. Comparison of compact scheme with analytical solution, sealed and killed 
end, point soma l:compact; :analytical 
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(c) N = 30 (d) TV = 40 



Figure 7. Comparison between compact, analytical and central difference schemes ( point 
soma), N = 10,20,30,40, ■:compact;H:analytical;B:central 



90 




(c) N = 30 (d) TV = 40 



Figure 8. Comparison of error between compact and central difference schemes(point 
soma), N = 10,20,30,40, ■icompact; :central 
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(a) sealed end (b) killed end 



Figure 9. Effect of diameter, sealed and killed end, (point soma), : d = 1.85* 10 4(cm); 

— : d = 3 * 1.85 * 1(T4; — -:d = 5 * 1.85 * 10~4; :d = 7 * 1.85 * 10"4; :d = 9 * 1.85 * 10"4 




(a) sealed end (b) killed end 



Figure 10. Evolution of voltage, current injected at i = /(point soma), sealed end and killed 
end boundaries, : n = 1; — : n = 100; :n == 500; :n == 2000; :n == 294670 
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(a) sealed end (b) killed end 



Figure 11. Current injected at i = N/2, (point soma) sealed end and killed end 
boundariesH:compact 
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Tables 

Table 1. Values of AT at various N 



N 


AT 




10 


8.8088e - 


04 


20 


1.9765e- 


04 


30 


8.4841e - 


05 


40 


4.6911e- 


05 


80 


1.1433e- 


05 



Table 2. Values of t in msec at various n 



n t 

1 070016968 

100 0.16968 

500 0.84840 

2000 3.3936 

294670 500 



Table 3. Parameters of dendrite used in simulation 



Parameter 


Values 


length 


400^ 


diameter 


2 * 1.85/x 




20 * I0 3 n.cm 2 




330S1. cm 


C m 


l^ifarad/cm 2 


T 


20 msec 


linj 


O.lnanoamperes 
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Table 4. Comparison of percentage error (Ep) between compact and central for N = 10(point soma) 



Error 


i = 1 


i = 2 


i = 3 


i = 4 


i = 5 


i = 6 


i = 7 


i = 8 


i = 9 


i = 10 


Compact 
Central 


0.098941 
0.097582 


0.096008 
0.094655 


0.093412 
0.092074 


0.091144 
0.089819 


0.089198 
0.087880 


0.087567 
0.086251 


0.086244 
0.084926 


0.085225 
0.083901 


0.084506 
0.083172 


0.084267 
0.082929 



Table 5. Resolving Efficiency e of the second derivative schemes, ( [8], Table 5) 



Scheme 


e = 0.1 


e = 0.01 


e = 0.001 


Fourth order central 


0.59 


0.31 


0.17 


Fourth order compact 


0.68 


0.39 


0.22 


Sixth order tridiagonal 


0.80 


0.55 


0.38 



Table 6. Values of A (cm) at various d(cm) 



d 




A 


1.85* 10-4 




0.052223 


3 * 1.85 * 10" 


4 


0.091701 


5 * 1.85 * 10" 


4 


0.11839 


7*1.85*10" 


4 


0.14008 


9 * 1.85 * 10" 


4 


0.15883 



